Setup

knitr::opts_chunk$set(echo = TRUE)

library(tidyverse)
library(R.utils)
library(wbCorr)
library(readxl)
library(kableExtra)
library(brms)
library(bayesplot)
library(beepr)
library(DHARMa)

source(file.path('Functions', 'ReportModels.R'))
source(file.path('Functions', 'PrettyTables.R'))
source(file.path('Functions', 'ReportMeasures.R'))
source(file.path('Functions', 'PrepareData.R'))
system("shutdown /a")
## [1] 1116
# Set options for analysis
use_mi = FALSE
shutdown = FALSE
report_ordinal = FALSE

options(
  dplyr.print_max = 100, 
  brms.backend = 'cmdstan',
  brms.file_refit = ifelse(use_mi, 'never', 'on_change'),
  brms.file_refit = 'always',
  error = function() beepr::beep(sound = 5)
)
df <- openxlsx::read.xlsx('long.xlsx')
df_original <- df

df_list <- prepare_data(df, use_mi = use_mi)

Constructing scales Re-coding pusing reshaping data (4field) centering data within and between

df_double <- df_list[[1]]
df_full <- df_list[[2]]

Descriptives

Sample Statistics

# Income (mean of both partner's report)
merge_income <- function(income1, income2) {
  merged_income <- numeric(length(income1))
  
  # Loop through each pair of incomes
  for (i in seq_along(income1)) {
    # Handle NA
    if (is.na(income1[i])) {
      merged_income[i] <- income2[i]
    } 
    
    else if (is.na(income2[i])) {
      merged_income[i] <- income1[i]
    }
    
    # if both are informative, take mean and round
    else if (income1[i] %in% 1:6 && income2[i] %in% 1:6) {
      merged_income[i] <- round((income1[i] + income2[i]) / 2)
    }
    
    # if one is informative and the other not, use the informative one. 
    else if (income1[i] %in% 1:6) {
      merged_income[i] <- income1[i]
    }
    else if (income2[i] %in% 1:6) {
      merged_income[i] <- income2[i]
    }
    
    # Now we only have cases, where both are either 70 or 99. We simply report "undisclosed". 
    else {
      merged_income[i] <- 99
    }
  }
  
  # Convert to factor
  merged_income <- factor(
    merged_income, levels = c(1,2,3,4,5,6,70,99), 
    labels = c(
      "up to CHF 2'000.-", 
      "CHF 2'001.- to CHF 4'000.-",
      "CHF 4'001.- to CHF 6'000.-",
      "CHF 6'001.- to CHF 8'000.-",
      "CHF 8'001.- to CHF 10'000.-", 
      "above CHF 10'000.-", 
      "I don't know",
      "Undisclosed"
  )
  )
  
  return(merged_income)
}



df_sample_report <- df_full %>%
  group_by(coupleID) %>%
  arrange(userID) %>%
  # Computing couple level variables
  mutate(
    Household_Income = merge_income(first(pre_income_1), last(pre_income_1)),
    
    reldur = pre_rel_duration_m / 12 + pre_rel_duration_y,
    Relationship_duration = mean(reldur, na.rm = TRUE),
    
    habdur = pre_hab_duration_m / 12 + pre_hab_duration_y,
    Cohabiting_duration = mean(habdur, na.rm = TRUE),
    
    Marital_status = factor(
      case_when(
        all(pre_mat_stat == 1) ~ "Married",
        any(pre_mat_stat == 1) ~ "One Partner Married",
        TRUE ~ "Not Married"
      )
    ),
    
    Have_children = factor(
      (first(pre_child_option) + last(pre_child_option)) > 0, 
      levels = c(FALSE, TRUE), 
      labels = c('Have Children', 'No Children')),
    
    Gender = factor(
      gender, 
      levels = c(1,2,3), 
      labels = c('Male','Female', 'Other')),
    

    Couple_type = as.factor(
      case_when(
        first(Gender) == last(Gender) & first(Gender) == 'Male' ~ 'Same-Sex Couple (Male)',
        first(Gender) == last(Gender) & first(Gender) == 'Female' ~ 'Same-Sex Couple (Female)',
        TRUE ~ 'Mixed-sex Couple'
      )
    )
  ) %>%
  ungroup() %>%
  # Individual level variables
  mutate(
    Age = pre_age,
    Handedness = factor(
      pre_handedness, 
      levels = c(0, 1, 2), 
      c('Right','Left', 'Ambidextrous')),
  Highest_Education = factor(
    pre_education, 
    levels = c(1,2,3,4,5,6,7), 
    labels = c(
      "(still) no school diploma",
      "compulsory education (9 years)",
      "vocational training (apprenticeship)",
      "Matura (university entrance qualification)",
      "Bachelor's degree", 
      "Master's degree",
      "Doctorate degree"
      )
    ),
  BMI = pre_weight / ((pre_height / 100)^2) # to meters
  ) %>%
  select(c(Relationship_duration, Cohabiting_duration, Couple_type, Household_Income, 
            Marital_status, Have_children, 
            Gender, Age, Handedness, Highest_Education, BMI))


sample_table <- report_measures(df_sample_report, ICC = F)
sample_table$n_Obs <- as.numeric(sample_table$n_Obs) / 55
rownames(sample_table) <- NULL

n_couple_vars <- 17
sample_table$n_Obs[1:n_couple_vars] <- sample_table$n_Obs[1:n_couple_vars] / 2

packing_sample <- list(
    "Couple level variables (38 couples)" = 
      c(1, n_couple_vars),
    "Individual level variables (76 individuals)" 
    = c(n_couple_vars+1, nrow(sample_table))
    ) 

df_sample_summary <- print_df(
  sample_table,
  rows_to_pack = packing_sample
)

export_xlsx(df_sample_summary,
            file.path('Output', 'SampleDescription.xlsx'),
            merge_option = 'both',
            rows_to_pack = packing_sample,
            colwidths = c(20,35,7,7,7,7,7,10)
            )
## 
## Attaching package: 'rvest'
## The following object is masked from 'package:readr':
## 
##     guess_encoding
df_sample_summary
Variable Level n_Obs percentage_Obs Missing Mean SD Range
Couple level variables (38 couples)
Relationship_duration NA 38 NA 0% 9.23 9.03 0.58-36.00
Cohabiting_duration NA 38 NA 0% 7.53 9.14 0.25-33.00
Couple_type Same-Sex Couple (Male) 1 3% NA NA NA NA
Couple_type Same-Sex Couple (Female) 1 3% NA NA NA NA
Couple_type Mixed-sex Couple 36 95% NA NA NA NA
Household_Income I don’t know 0 0% NA NA NA NA
Household_Income up to CHF 2’000.- 2 5% NA NA NA NA
Household_Income CHF 2’001.- to CHF 4’000.- 3 8% NA NA NA NA
Household_Income CHF 6’001.- to CHF 8’000.- 3 8% NA NA NA NA
Household_Income CHF 4’001.- to CHF 6’000.- 5 13% NA NA NA NA
Household_Income Undisclosed 5 13% NA NA NA NA
Household_Income CHF 8’001.- to CHF 10’000.- 8 21% NA NA NA NA
Household_Income above CHF 10’000.- 12 32% NA NA NA NA
Marital_status Married 13 34% NA NA NA NA
Marital_status Not Married 25 66% NA NA NA NA
Have_children No Children 11 29% NA NA NA NA
Have_children Have Children 27 71% NA NA NA NA
Individual level variables (76 individuals)
Gender Other 0 0% NA NA NA NA
Gender Male 38 50% NA NA NA NA
Gender Female 38 50% NA NA NA NA
Age NA 76 NA 0% 34.01 10.96 19.00-60.00
Handedness Ambidextrous 1 1% NA NA NA NA
Handedness Left 11 14% NA NA NA NA
Handedness Right 64 84% NA NA NA NA
Highest_Education (still) no school diploma 0 0% NA NA NA NA
Highest_Education compulsory education (9 years) 0 0% NA NA NA NA
Highest_Education Doctorate degree 0 0% NA NA NA NA
Highest_Education vocational training (apprenticeship) 8 11% NA NA NA NA
Highest_Education Master’s degree 21 28% NA NA NA NA
Highest_Education Bachelor’s degree 23 30% NA NA NA NA
Highest_Education Matura (university entrance qualification) 24 32% NA NA NA NA
BMI NA 76 NA 0% 24.94 4.08 16.37-33.95

Measures

Focal Variables

main_constructs <- c("persuasion", "pressure","pushing", 
                     "pa_sub", "pa_obj", "aff", "reactance"
                     )

main_descriptives <- report_measures(
  data = df_full, 
  measures = main_constructs,
  ICC = TRUE, 
  cluster_var = df_full$userID)

openxlsx::write.xlsx(
  main_descriptives, 
  file.path('Output', 'DescriptivesMain.xlsx')
  )

print_df(main_descriptives)
Variable n_Obs Missing Mean SD Range ICC
persuasion 4180 6% 0.42 1.08 0.00- 5.00 0.23
pressure 4180 6% 0.12 0.62 0.00- 5.00 0.58
pushing 4180 6% 0.16 0.66 0.00- 5.00 0.11
pa_sub 4180 6% 30.40 54.78 0.00-720.00 0.16
pa_obj 4180 11% 144.41 117.81 5.75-971.25 0.18
aff 4180 6% 4.83 1.14 1.00- 6.00 0.41
reactance 4180 81% 0.79 1.32 0.00- 5.00 0.47

All Variables

all_constructs <- c(
  main_constructs,
  "day",
  "weartime",
  "isWeekend",
  "plan",
  "studyGroup",
  "support",
  "got_JITAI",
  "skilled_support"
)


all_descriptives <- report_measures(df_full, all_constructs, ICC = F)

openxlsx::write.xlsx(
  all_descriptives, 
  file.path('Output', 'DescriptivesAll.xlsx')
)

print_df(all_descriptives)
Variable Level n_Obs percentage_Obs Missing Mean SD Range
persuasion NA 4180 NA 6% 0.42 1.08 0.00- 5.00
pressure NA 4180 NA 6% 0.12 0.62 0.00- 5.00
pushing NA 4180 NA 6% 0.16 0.66 0.00- 5.00
pa_sub NA 4180 NA 6% 30.40 54.78 0.00- 720.00
pa_obj NA 4180 NA 11% 144.41 117.81 5.75- 971.25
aff NA 4180 NA 6% 4.83 1.14 1.00- 6.00
reactance NA 4180 NA 81% 0.79 1.32 0.00- 5.00
day NA 4180 NA 0% 0.50 0.29 0.00- 1.00
weartime NA 4180 NA 1% 1057.42 384.29 0.00-1440.00
isWeekend Weekend 1216 29% NA NA NA NA
isWeekend Weekday 2964 71% NA NA NA NA
plan missing 238 6% NA NA NA NA
plan Plan 1860 44% NA NA NA NA
plan No plan 2082 50% NA NA NA NA
studyGroup last 3 weeks interventions 1320 32% NA NA NA NA
studyGroup Allways inerventions 1430 34% NA NA NA NA
studyGroup First 3 weeks interventions 1430 34% NA NA NA NA
support NA 4180 NA 6% 0.91 1.52 0.00- 5.00
got_JITAI JITAI received 585 14% NA NA NA NA
got_JITAI No JITAI 3595 86% NA NA NA NA
skilled_support Days before Intervention 1036 25% NA NA NA NA
skilled_support Days after Intervention 3144 75% NA NA NA NA

Correlations

cors <- wbCorr(df_full[,c(main_constructs)], df_full$coupleID, method = 'spearman')

main_cors <- summary(cors, 'wb')$merged_wb


openxlsx::write.xlsx(
  main_cors, 
  file.path('Output', 'Correlations.xlsx')
)

print_df(main_cors, width = '7em')
persuasion pressure pushing pa_sub pa_obj aff reactance
persuasion [0.20] 0.21*** 0.40*** 0.17*** 0.07*** 0.00 -0.05
pressure 0.30 [0.55] 0.28*** -0.03 -0.04* -0.01 0.26***
pushing 0.63*** 0.40* [0.08] 0.14*** 0.05** 0.01 0.07*
pa_sub 0.27 -0.10 0.24 [0.15] 0.31*** 0.20*** -0.04
pa_obj 0.13 -0.08 0.27 0.51*** [0.13] 0.19*** 0.06
aff 0.28 -0.07 0.29 0.52*** 0.22 [0.28] -0.01
reactance 0.18 0.23 0.11 0.06 0.38* -0.12 [0.42]

Within-person correlations are above the diagonal and between-person correlations are below the diagonal.
On the diagonal are intraclass correlations (ICCs)

Modelling

# For indistinguishable Dyads
model_rows_fixed <- c(
    'Intercept', 
    # '-- WITHIN PERSON MAIN EFFECTS --', 
    'persuasion_self_cw', 
    'persuasion_partner_cw', 
    'pressure_self_cw', 
    'pressure_partner_cw', 
    'pushing_self_cw', 
    'pushing_partner_cw', 
    'day', 
    'weartime_self_cw',
    
    # '-- BETWEEN PERSON MAIN EFFECTS',
    'persuasion_self_cb',
    'persuasion_partner_cb',
    'pressure_self_cb',
    'pressure_partner_cb',
    'pushing_self_cb',
    'pushing_partner_cb',
    'weartime_self_cb'
  )


model_rows_fixed_ordinal <- c(
  model_rows_fixed[1],
  'Intercept[1]',
  'Intercept[2]',
  'Intercept[3]',
  'Intercept[4]',
  'Intercept[5]',
  model_rows_fixed[2:length(model_rows_fixed)]
)

model_rows_random <- c(
  # '--------------',
  # '-- RANDOM EFFECTS --',
  'sd(Intercept)', 
  'sd(persuasion_self_cw)',
  'sd(persuasion_partner_cw)',
  'sd(pressure_self_cw)',
  'sd(pressure_partner_cw)',
  'sd(pushing_self_cw)',
  'sd(pushing_partner_cw)',
  # '-- CORRELATION STRUCTURE -- ', 
  'ar[1]', 
  'nu',
  'shape',
  'sderr',
  'sigma'
)

model_rows_random_ordinal <- c(model_rows_random,'disc')
# For indistinguishable Dyads
model_rownames_fixed <- c(
    'Intercept', 
    # '-- WITHIN PERSON MAIN EFFECTS --', 
    'Daily received persuasion target -> target', 
    'Daily received persuasion target -> agent', 
    'Daily received pressure target -> target', 
    'Daily received pressure target -> agent', 
    'Daily received pushing target -> target', 
    'Daily received pushing target -> agent', 
    'Day', 
    'Daily weartime',
    
    # '-- BETWEEN PERSON MAIN EFFECTS',
    'Mean received persuasion target -> target',
    'Mean received persuasion target -> agent',
    'Mean received pressure target -> target',
    'Mean received pressure target -> agent',
    'Mean received pushing target -> target',
    'Mean received pushing target -> agent',
    'Mean weartime'
  )


model_rownames_fixed_ordinal <- c(
  model_rownames_fixed[1],
  'Intercept[1]',
  'Intercept[2]',
  'Intercept[3]',
  'Intercept[4]',
  'Intercept[5]',
  model_rownames_fixed[2:length(model_rownames_fixed)]
)

model_rownames_random <- c(
  # '--------------',
  # '-- RANDOM EFFECTS --',
  'sd(Intercept)', 
  'sd(Daily received persuasion target -> target)', 
  'sd(Daily received persuasion target -> agent)', 
  'sd(Daily received pressure target -> target)', 
  'sd(Daily received pressure target -> agent)', 
  'sd(Daily received pushing target -> target)', 
  'sd(Daily received pushing target -> agent)', 
  # '-- CORRELATION STRUCTURE -- ', 
  'ar[1]', 
  'nu',
  'shape',
  'sderr',
  'sigma'
)

model_rownames_random_ordinal <- c(model_rownames_random,'disc')
# For indistinguishable Dyads
model_rownames_fixed <- c(
    "Intercept", 
    # "-- WITHIN PERSON MAIN EFFECTS --", 
    "Daily persuasion experienced", 
    "Daily persuasion utilized (partner's view)", # OR partner received
    "Daily pressure experienced", 
    "Daily pressure utilized (partner's view)", 
    "Daily pushing experienced", 
    "Daily pushing utilized (partner's view)", 
    "Day", 
    "Daily weartime",
    
    # "-- BETWEEN PERSON MAIN EFFECTS",
    "Mean persuasion experienced", 
    "Mean persuasion utilized (partner's view)", 
    "Mean pressure experienced", 
    "Mean pressure utilized (partner's view)", 
    "Mean pushing experienced", 
    "Mean pushing utilized (partner's view)", 
    "Mean weartime"
  )


model_rownames_fixed_ordinal <- c(
  model_rownames_fixed[1],
  'Intercept[1]',
  'Intercept[2]',
  'Intercept[3]',
  'Intercept[4]',
  'Intercept[5]',
  model_rownames_fixed[2:length(model_rownames_fixed)]
)

model_rownames_random <- c(
  # '--------------',
  # '-- RANDOM EFFECTS --',
  'sd(Intercept)', 
  "sd(Daily persuasion experienced)", 
  "sd(Daily persuasion utilized (partner's view))", # OR partner received
  "sd(Daily pressure experienced)", 
  "sd(Daily pressure utilized (partner's view))", 
  "sd(Daily pushing experienced)", 
  "sd(Daily pushing utilized (partner's view))", 
  # '-- CORRELATION STRUCTURE -- ', 
  'ar[1]', 
  'nu',
  'shape',
  'sderr',
  'sigma'
)

model_rownames_random_ordinal <- c(model_rownames_random,'disc')
rows_to_pack <- list(
  "Within-Person Effects" = c(2,9),
  "Between-Person Effects" = c(10,16),
  "Random Effects" = c(17, 23), 
  "Additional Parameters" = c(24,28)
  )


rows_to_pack_ordinal <- list(
  "Intercepts" = c(1,6),
  "Within-Person Effects" = c(2+5,9+5),
  "Between-Person Effects" = c(10+5,16+5),
  "Random Effects" = c(17+5, 23+5), 
  "Additional Parameters" = c(24+5,28+6)
  )

Subjective MVPA

range(df_double$pa_sub, na.rm = T) 
## [1]   0 720
hist(df_double$pa_sub, breaks = 100)

Modelling using the gaussian family fails. Due to the many zeros, transformations won’t help estimating the models. We employ the negative binomial family.

formula <- bf(
  pa_sub ~ 
    persuasion_self_cw + persuasion_partner_cw +
    pressure_self_cw + pressure_partner_cw +
    pushing_self_cw + pushing_partner_cw +
    
    persuasion_self_cb + persuasion_partner_cb +
    pressure_self_cb + pressure_partner_cb +
    pushing_self_cb + pushing_partner_cb +
    
    day + 
    
    # Random effects
    (persuasion_self_cw + persuasion_partner_cw +
    pressure_self_cw + pressure_partner_cw +
    pushing_self_cw + pushing_partner_cw | coupleID),

  autocor = ~ ar(time = day, gr = coupleID:userID, p = 1)
)




prior1 <- c(
  brms::set_prior("normal(0, 2.5)", class = "b"),
  brms::set_prior("normal(0, 50)", class = "Intercept", lb = 0),
  
  brms::set_prior("normal(0, 10)", class = "sd", group = "coupleID", lb = 0),
  
  brms::set_prior("cauchy(0, 5)", class = "ar", lb = -1, ub = 1),
  brms::set_prior("cauchy(0, 20)", class = "shape"), 
  brms::set_prior("cauchy(0, 10)", class='sderr')
)


#df_minimal <- df_double[, c("AorB", all.vars(as.formula(formula)))]

pa_sub <- my_brm(
  mi = use_mi, 
  imputed_data = implist,
  
  formula = formula, 
  prior = prior1,
  data = df_double, 
  family = brms::negbinomial(),
  #control = list(adapt_delta = 0.99),
  iter = 5000,
  warmup = 2000,
  chains = 4,
  cores = 4,
  seed = 7777,
  file = file.path("models_cache_brms", "pa_sub")
)
plot(pa_sub, ask = FALSE)

plot(pp_check(pa_sub, type = 'ecdf_overlay'))
## Using 10 posterior draws for ppc type 'ecdf_overlay' by default.

plot(pp_check(pa_sub))
## Using 10 posterior draws for ppc type 'dens_overlay' by default.

pp_check_transformed(pa_sub, transform = log1p)

loo(pa_sub)
## 
## Computed from 12000 by 3736 log-likelihood matrix.
## 
##          Estimate    SE
## elpd_loo -12068.5 177.4
## p_loo        37.8   3.0
## looic     24136.9 354.9
## ------
## MCSE of elpd_loo is NA.
## MCSE and ESS estimates assume MCMC draws (r_eff in [0.3, 2.0]).
## 
## Pareto k diagnostic values:
##                          Count Pct.    Min. ESS
## (-Inf, 0.7]   (good)     3730  99.8%   1241    
##    (0.7, 1]   (bad)         6   0.2%   <NA>    
##    (1, Inf)   (very bad)    0   0.0%   <NA>    
## See help('pareto-k-diagnostic') for details.
summarize_brms(
  pa_sub, 
  model_rows_fixed = model_rows_fixed,
  model_rows_random = model_rows_random,
  model_rownames_fixed = model_rownames_fixed,
  model_rownames_random = model_rownames_random,
  exponentiate = T) %>%
  print_df(rows_to_pack = rows_to_pack)
IRR l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
Intercept 24.25* 17.04 34.72 1.002 1749.28 3624.21
Within-Person Effects
Daily persuasion experienced 0.97 0.87 1.07 1.000 14143.64 9992.62
Daily persuasion utilized (partner’s view) 1.05 0.95 1.17 1.001 13107.42 9379.26
Daily pressure experienced 1.12 0.89 1.45 1.001 13954.29 9470.29
Daily pressure utilized (partner’s view) 0.99 0.79 1.28 1.000 13827.21 8061.82
Daily pushing experienced 0.93 0.80 1.09 1.001 13565.74 8610.64
Daily pushing utilized (partner’s view) 0.92 0.79 1.09 1.000 14123.42 9220.05
Day 1.30 0.95 1.78 1.000 16336.30 9775.92
Daily weartime NA NA NA NA NA NA
Between-Person Effects
Mean persuasion experienced 1.21 0.72 2.06 1.000 5760.83 7560.12
Mean persuasion utilized (partner’s view) 1.16 0.71 1.94 1.000 5629.19 8117.99
Mean pressure experienced 1.60 0.82 3.12 1.000 6742.86 8277.70
Mean pressure utilized (partner’s view) 0.53 0.28 1.04 1.000 8351.55 9243.99
Mean pushing experienced 0.43* 0.19 0.99 1.000 7157.94 8917.45
Mean pushing utilized (partner’s view) 0.68 0.29 1.55 1.000 7035.02 9067.82
Mean weartime NA NA NA NA NA NA
Random Effects
sd(Intercept) 0.68 0.51 0.91 1.00 3296.67 6251.21
sd(Daily persuasion experienced) 0.22 0.07 0.37 1.00 3517.06 2486.89
sd(Daily persuasion utilized (partner’s view)) 0.18 0.03 0.33 1.00 3661.30 2705.55
sd(Daily pressure experienced) 0.14 0.01 0.44 1.00 7508.04 5234.97
sd(Daily pressure utilized (partner’s view)) 0.15 0.01 0.45 1.00 7110.00 6357.41
sd(Daily pushing experienced) 0.25 0.02 0.53 1.00 3541.13 3656.25
sd(Daily pushing utilized (partner’s view)) 0.14 0.01 0.36 1.00 4790.81 4247.33
Additional Parameters
ar[1] 0.02 -0.93 0.94 1.00 12676.91 7882.74
nu NA NA NA NA NA NA
shape 0.14 0.13 0.14 1.00 15982.23 9076.33
sderr 0.05 0.00 0.13 1.00 6828.15 5892.48
sigma NA NA NA NA NA NA

Device Based MVPA

range(df_double$pa_obj, na.rm = T) 
## [1]   5.75 971.25
hist(df_double$pa_obj, breaks = 50)

df_double$pa_obj_log <- log(df_double$pa_obj)

hist(df_double$pa_obj_log, breaks = 50)

We tried negative binomial here as well for consistency, but the model did not converge. Poisson also did not work. As we have no zeros in this distribution, we log transform.

formula <- bf(
  pa_obj_log ~ 
    persuasion_self_cw + persuasion_partner_cw +
    pressure_self_cw + pressure_partner_cw +
    pushing_self_cw + pushing_partner_cw +
    
    persuasion_self_cb + persuasion_partner_cb +
    pressure_self_cb + pressure_partner_cb +
    pushing_self_cb + pushing_partner_cb +
    
    day + weartime_self_cw + weartime_self_cb +
    
    # Random effects
    (persuasion_self_cw + persuasion_partner_cw +
    pressure_self_cw + pressure_partner_cw +
    pushing_self_cw + pushing_partner_cw | coupleID),

  autocor = ~ ar(time = day, gr = coupleID:userID, p = 1)
)



prior1 <- c(
  brms::set_prior("normal(0, 2.5)", class = "b"),
  brms::set_prior("normal(0, 50)", class = "Intercept", lb = 0),
  
  brms::set_prior("normal(0, 10)", class = "sd", group = "coupleID", lb = 0),
  
  brms::set_prior("cauchy(0, 5)", class = "ar", lb = -1, ub = 1),
  brms::set_prior("cauchy(0, 10)", class = "sigma", lb = 0)
)


#df_minimal <- df_double[, c("AorB", all.vars(as.formula(formula)))]

pa_obj_log <- my_brm(
  mi = use_mi, 
  imputed_data = implist,
  
  formula = formula, 
  prior = prior1,
  data = df_double, 
  family = gaussian(),
  #control = list(adapt_delta = 0.99),
  iter = 5000,
  warmup = 2000,
  chains = 4,
  cores = 4,
  seed = 7777,
  file = file.path("models_cache_brms", "pa_obj_log")
)
plot(pa_obj_log, ask = FALSE)

plot(pp_check(pa_obj_log, type = 'ecdf_overlay'))
## Using 10 posterior draws for ppc type 'ecdf_overlay' by default.

plot(pp_check(pa_obj_log))
## Using 10 posterior draws for ppc type 'dens_overlay' by default.

#pp_check_transformed(pa_obj_log, transform = log1p)

loo(pa_obj_log)
## 
## Computed from 12000 by 3337 log-likelihood matrix.
## 
##          Estimate    SE
## elpd_loo  -2823.4  55.5
## p_loo       101.8   4.8
## looic      5646.8 111.1
## ------
## MCSE of elpd_loo is NA.
## MCSE and ESS estimates assume MCMC draws (r_eff in [0.3, 2.4]).
## 
## Pareto k diagnostic values:
##                          Count Pct.    Min. ESS
## (-Inf, 0.7]   (good)     3336  100.0%  916     
##    (0.7, 1]   (bad)         1    0.0%  <NA>    
##    (1, Inf)   (very bad)    0    0.0%  <NA>    
## See help('pareto-k-diagnostic') for details.
summarize_brms(
  pa_obj_log, 
  model_rows_fixed = model_rows_fixed,
  model_rows_random = model_rows_random,
  model_rownames_fixed = model_rownames_fixed,
  model_rownames_random = model_rownames_random,
  exponentiate = T) %>%
  print_df(rows_to_pack = rows_to_pack)
exp(Est.) l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
Intercept 116.98* 104.97 130.74 1.000 2134.78 4489.10
Within-Person Effects
Daily persuasion experienced 1.00 0.98 1.02 1.000 20457.76 9005.79
Daily persuasion utilized (partner’s view) 1.00 0.98 1.02 1.000 18658.71 8945.50
Daily pressure experienced 1.03 0.98 1.09 1.000 18727.08 9716.36
Daily pressure utilized (partner’s view) 0.98 0.93 1.03 1.000 21135.85 9438.79
Daily pushing experienced 1.00 0.97 1.03 1.000 21479.39 9082.80
Daily pushing utilized (partner’s view) 1.02 0.99 1.06 1.001 18756.07 9694.70
Day 1.03 0.94 1.12 1.001 21218.74 8889.38
Daily weartime 1.00 1.00 1.00 1.001 12170.05 8378.04
Between-Person Effects
Mean persuasion experienced 1.05 0.90 1.22 1.000 4211.44 6575.30
Mean persuasion utilized (partner’s view) 1.05 0.90 1.23 1.000 4152.95 7041.28
Mean pressure experienced 0.86 0.72 1.02 1.000 7520.72 8624.58
Mean pressure utilized (partner’s view) 0.87 0.74 1.01 1.000 6678.31 8130.53
Mean pushing experienced 1.12 0.90 1.40 1.001 6417.14 8251.81
Mean pushing utilized (partner’s view) 1.12 0.90 1.38 1.002 6206.59 7988.07
Mean weartime 1.00 1.00 1.00 1.000 12789.09 9738.02
Random Effects
sd(Intercept) 0.29 0.22 0.38 1.00 2763.57 4596.67
sd(Daily persuasion experienced) 0.06 0.03 0.09 1.00 8017.53 8058.29
sd(Daily persuasion utilized (partner’s view)) 0.05 0.02 0.08 1.00 6347.96 5412.62
sd(Daily pressure experienced) 0.05 0.00 0.14 1.00 4468.00 5361.23
sd(Daily pressure utilized (partner’s view)) 0.04 0.00 0.11 1.00 7259.06 5861.92
sd(Daily pushing experienced) 0.08 0.01 0.15 1.00 2846.80 2766.55
sd(Daily pushing utilized (partner’s view)) 0.04 0.00 0.10 1.00 3804.63 5277.03
Additional Parameters
ar[1] 0.28 0.25 0.32 1.00 20704.78 8314.00
nu NA NA NA NA NA NA
shape NA NA NA NA NA NA
sderr NA NA NA NA NA NA
sigma 0.55 0.54 0.57 1.00 20979.44 8865.71

Affect

range(df_double$aff, na.rm = T) 
## [1] 1 6
hist(df_double$aff, breaks = 15)

formula <- bf(
  aff ~ 
    persuasion_self_cw + persuasion_partner_cw +
    pressure_self_cw + pressure_partner_cw +
    pushing_self_cw + pushing_partner_cw +
    
    persuasion_self_cb + persuasion_partner_cb +
    pressure_self_cb + pressure_partner_cb +
    pushing_self_cb + pushing_partner_cb +
    
    day +
    
    # Random effects
    (persuasion_self_cw + persuasion_partner_cw +
    pressure_self_cw + pressure_partner_cw +
    pushing_self_cw + pushing_partner_cw | coupleID),

  autocor = ~ ar(time = day, gr = coupleID:userID, p = 1)
)



prior1 <- c(
  brms::set_prior("normal(0, 5)", class = "b"),
  brms::set_prior("normal(0, 20)", class = "Intercept", lb=0, ub=6), # range of the outcome scale
  
  brms::set_prior("normal(0, 2)", class = "sd", group = "coupleID", lb = 0),

  brms::set_prior("cauchy(0, 5)", class = "ar", lb = -1, ub = 1),
  brms::set_prior("cauchy(0, 10)", class = "sigma", lb = 0)
  
)

#df_minimal <- df_double[, c("AorB", all.vars(as.formula(formula)))]

mood <- my_brm(
  mi = use_mi, 
  imputed_data = implist,
  
  formula = formula, 
  prior = prior1,
  data = df_double, 
  family = gaussian(),
  #control = list(adapt_delta = 0.95, max_treedepth = 15),
  iter = 5000,
  warmup = 2000,
  chains = 4,
  cores = 4,
  seed = 7777,
  file = file.path("models_cache_brms", "mood")
)
plot(mood, ask = FALSE)

plot(pp_check(mood, type = 'ecdf_overlay'))
## Using 10 posterior draws for ppc type 'ecdf_overlay' by default.

plot(pp_check(mood))
## Using 10 posterior draws for ppc type 'dens_overlay' by default.

#pp_check_transformed(mood, transform = log1p)

loo(pa_sub)
## 
## Computed from 12000 by 3736 log-likelihood matrix.
## 
##          Estimate    SE
## elpd_loo -12068.5 177.4
## p_loo        37.8   3.0
## looic     24136.9 354.9
## ------
## MCSE of elpd_loo is NA.
## MCSE and ESS estimates assume MCMC draws (r_eff in [0.3, 2.0]).
## 
## Pareto k diagnostic values:
##                          Count Pct.    Min. ESS
## (-Inf, 0.7]   (good)     3730  99.8%   1241    
##    (0.7, 1]   (bad)         6   0.2%   <NA>    
##    (1, Inf)   (very bad)    0   0.0%   <NA>    
## See help('pareto-k-diagnostic') for details.
summarize_brms(
  mood, 
  model_rows_fixed = model_rows_fixed,
  model_rows_random = model_rows_random,
  model_rownames_fixed = model_rownames_fixed,
  model_rownames_random = model_rownames_random,
  exponentiate = F) %>%
  print_df(rows_to_pack = rows_to_pack)
b l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
Intercept 4.93* 4.71 5.15 1.002 1729.60 4045.66
Within-Person Effects
Daily persuasion experienced 0.00 -0.03 0.03 1.000 17062.16 9995.68
Daily persuasion utilized (partner’s view) -0.01 -0.05 0.02 1.000 16978.58 9218.48
Daily pressure experienced 0.00 -0.08 0.07 1.000 16091.41 9515.63
Daily pressure utilized (partner’s view) 0.05 -0.03 0.13 1.000 17293.06 9336.91
Daily pushing experienced -0.03 -0.08 0.02 1.000 16952.08 9530.67
Daily pushing utilized (partner’s view) -0.04 -0.09 0.00 1.000 17521.29 9489.19
Day -0.10 -0.25 0.06 1.001 21634.92 9437.78
Daily weartime NA NA NA NA NA NA
Between-Person Effects
Mean persuasion experienced 0.14 -0.13 0.42 1.001 3845.63 6120.15
Mean persuasion utilized (partner’s view) 0.21 -0.06 0.49 1.001 3728.88 5373.73
Mean pressure experienced -0.11 -0.37 0.14 1.000 5311.75 7817.81
Mean pressure utilized (partner’s view) -0.02 -0.28 0.23 1.000 5214.22 7073.42
Mean pushing experienced -0.20 -0.60 0.21 1.000 5631.69 7635.04
Mean pushing utilized (partner’s view) -0.43* -0.83 -0.03 1.000 5605.88 7625.27
Mean weartime NA NA NA NA NA NA
Random Effects
sd(Intercept) 0.63 0.49 0.81 1.00 2995.72 5082.68
sd(Daily persuasion experienced) 0.03 0.00 0.07 1.00 4868.91 5412.03
sd(Daily persuasion utilized (partner’s view)) 0.06 0.01 0.11 1.00 2844.71 3592.91
sd(Daily pressure experienced) 0.09 0.00 0.26 1.00 4479.26 5751.81
sd(Daily pressure utilized (partner’s view)) 0.15 0.01 0.33 1.00 3639.77 4174.77
sd(Daily pushing experienced) 0.09 0.02 0.16 1.00 5226.28 3542.48
sd(Daily pushing utilized (partner’s view)) 0.08 0.01 0.16 1.00 4014.82 4191.50
Additional Parameters
ar[1] 0.45 0.42 0.48 1.00 19466.05 9345.83
nu NA NA NA NA NA NA
shape NA NA NA NA NA NA
sderr NA NA NA NA NA NA
sigma 0.87 0.85 0.89 1.00 20700.26 8273.28

Reactance

Gaussian

range(df_double$reactance, na.rm = T) 
## [1] 0 5
hist(df_double$reactance, breaks = 6)

formula <- bf(
  reactance ~ 
    persuasion_self_cw + persuasion_partner_cw +
    pressure_self_cw + pressure_partner_cw +
    pushing_self_cw + pushing_partner_cw +
    
    persuasion_self_cb + persuasion_partner_cb +
    pressure_self_cb + pressure_partner_cb +
    pushing_self_cb + pushing_partner_cb +
    
    day +
    
    # Random effects
    (persuasion_self_cw + persuasion_partner_cw +
    pressure_self_cw + pressure_partner_cw +
    pushing_self_cw + pushing_partner_cw | coupleID),

  autocor = ~ ar(time = day, gr = coupleID:userID, p = 1)
)



prior1 <- c(
  brms::set_prior("normal(0, 5)", class = "b"),
  brms::set_prior("normal(0, 20)", class = "Intercept", lb=0, ub=5), # range of the outcome scale
  
  brms::set_prior("normal(0, 2)", class = "sd", group = "coupleID", lb = 0),

  brms::set_prior("cauchy(0, 5)", class = "ar", lb = -1, ub = 1),
  brms::set_prior("cauchy(0, 10)", class = "sigma", lb = 0)
)




#df_minimal <- df_double[, c("AorB", all.vars(as.formula(formula)))]

reactance <- my_brm(
  mi = use_mi, 
  imputed_data = implist,
  
  formula = formula, 
  prior = prior1,
  data = df_double, 
  family = gaussian(),
  #control = list(adapt_delta = 0.95, max_treedepth = 15),
  iter = 5000,
  warmup = 2000,
  chains = 4,
  cores = 4,
  seed = 7777,
  file = file.path("models_cache_brms", "reactance")
)
plot(reactance, ask = FALSE)

plot(pp_check(reactance, type = 'ecdf_overlay'))
## Using 10 posterior draws for ppc type 'ecdf_overlay' by default.

plot(pp_check(reactance))
## Using 10 posterior draws for ppc type 'dens_overlay' by default.

#pp_check_transformed(reactance, transform = log1p)

loo(reactance)
## 
## Computed from 12000 by 756 log-likelihood matrix.
## 
##          Estimate   SE
## elpd_loo  -1074.7 34.0
## p_loo        85.8  8.2
## looic      2149.3 68.1
## ------
## MCSE of elpd_loo is NA.
## MCSE and ESS estimates assume MCMC draws (r_eff in [0.5, 1.7]).
## 
## Pareto k diagnostic values:
##                          Count Pct.    Min. ESS
## (-Inf, 0.7]   (good)     747   98.8%   208     
##    (0.7, 1]   (bad)        9    1.2%   <NA>    
##    (1, Inf)   (very bad)   0    0.0%   <NA>    
## See help('pareto-k-diagnostic') for details.
summarize_brms(
  reactance, 
  model_rows_fixed = model_rows_fixed,
  model_rows_random = model_rows_random,
  model_rownames_fixed = model_rownames_fixed,
  model_rownames_random = model_rownames_random,
  exponentiate = F) %>%
  print_df(rows_to_pack = rows_to_pack)
b l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
Intercept 0.64* 0.37 0.90 1.000 4289.28 5359.99
Within-Person Effects
Daily persuasion experienced 0.01 -0.04 0.06 1.000 17210.46 9290.26
Daily persuasion utilized (partner’s view) -0.04 -0.11 0.02 1.001 18013.65 9361.93
Daily pressure experienced -0.04 -0.13 0.06 1.000 19478.42 9281.07
Daily pressure utilized (partner’s view) -0.04 -0.17 0.09 1.000 20016.91 8704.37
Daily pushing experienced -0.01 -0.07 0.05 1.000 19261.80 8591.46
Daily pushing utilized (partner’s view) 0.02 -0.06 0.10 1.000 18358.77 9275.06
Day -0.25 -0.52 0.03 1.000 17817.26 10160.99
Daily weartime NA NA NA NA NA NA
Between-Person Effects
Mean persuasion experienced -0.04 -0.39 0.31 1.000 6477.01 8903.37
Mean persuasion utilized (partner’s view) 0.09 -0.30 0.47 1.000 7524.57 8543.92
Mean pressure experienced 0.29 -0.07 0.65 1.000 10033.44 9278.46
Mean pressure utilized (partner’s view) -0.22 -0.60 0.16 1.000 9827.60 8802.26
Mean pushing experienced -0.24 -0.77 0.30 1.000 8717.20 8554.84
Mean pushing utilized (partner’s view) 0.03 -0.56 0.60 1.000 8437.21 9210.15
Mean weartime NA NA NA NA NA NA
Random Effects
sd(Intercept) 0.57 0.43 0.76 1.00 5575.08 7441.66
sd(Daily persuasion experienced) 0.06 0.00 0.14 1.00 3401.29 6310.48
sd(Daily persuasion utilized (partner’s view)) 0.04 0.00 0.11 1.00 5881.87 4964.89
sd(Daily pressure experienced) 0.44 0.28 0.67 1.00 7323.82 8035.68
sd(Daily pressure utilized (partner’s view)) 0.22 0.02 0.58 1.00 3224.98 4830.01
sd(Daily pushing experienced) 0.09 0.00 0.23 1.00 2342.57 4591.05
sd(Daily pushing utilized (partner’s view)) 0.04 0.00 0.14 1.00 7840.01 8176.19
Additional Parameters
ar[1] 0.01 -0.07 0.10 1.00 13020.01 9177.59
nu NA NA NA NA NA NA
shape NA NA NA NA NA NA
sderr NA NA NA NA NA NA
sigma 0.94 0.89 0.99 1.00 12075.95 8745.61 
hypothesis(reactance, "pressure_self_cw > pushing_self_cw")
## Hypothesis Tests for class b:
##                 Hypothesis Estimate Est.Error CI.Lower CI.Upper Evid.Ratio
## 1 (pressure_self_cw... > 0    -0.03      0.06    -0.13     0.07       0.51
##   Post.Prob Star
## 1      0.34     
## ---
## 'CI': 90%-CI for one-sided and 95%-CI for two-sided hypotheses.
## '*': For one-sided hypotheses, the posterior probability exceeds 95%;
## for two-sided hypotheses, the value tested against lies outside the 95%-CI.
## Posterior probabilities of point hypotheses assume equal prior probabilities.

Ordinal

df_double$reactance_ordinal <- factor(df_double$reactance,
                                      levels = 0:5, 
                                      ordered = TRUE)

formula <- bf(
  reactance_ordinal ~ 
    persuasion_self_cw + persuasion_partner_cw +
    pressure_self_cw + pressure_partner_cw +
    pushing_self_cw + pushing_partner_cw +
    
    persuasion_self_cb + persuasion_partner_cb +
    pressure_self_cb + pressure_partner_cb +
    pushing_self_cb + pushing_partner_cb +
    
    day +
    
    # Random effects
    (persuasion_self_cw + persuasion_partner_cw +
    pressure_self_cw + pressure_partner_cw +
    pushing_self_cw + pushing_partner_cw | coupleID),

  autocor = ~ ar(time = day, gr = coupleID:userID, p = 1)
)



prior1 <- c(
  set_prior("normal(0, 2.5)", class = "b"),
  set_prior("normal(0, 5)", class = "Intercept"),
  set_prior("normal(0, 1)", class = "sd", group = "coupleID", lb = 0),
  set_prior("normal(0, 0.075)", class = "ar", lb = -1, ub = 1),
  set_prior("normal(0.5, 2.0)", class = "sderr", lb = 0)
)




#df_minimal <- df_double[, c("AorB", all.vars(as.formula(formula)))]

reactance_ordinal <- my_brm(
  mi = use_mi, 
  imputed_data = implist,
  
  formula = formula, 
  prior = prior1,
  data = df_double, 
  family = brms::cumulative(),
  control = list(adapt_delta = 0.95),
  iter = 5000,
  warmup = 2000,
  chains = 4,
  cores = 4,
  seed = 7777,
  file = file.path("models_cache_brms", "reactance_ordinal")
)
plot(reactance_ordinal, ask = FALSE)

plot(pp_check(reactance_ordinal, type = 'ecdf_overlay'))
## Using 10 posterior draws for ppc type 'ecdf_overlay' by default.

plot(pp_check(reactance_ordinal))
## Using 10 posterior draws for ppc type 'dens_overlay' by default.

#pp_check_transformed(reactance_ordinal, transform = log1p)

loo(reactance_ordinal)
## 
## Computed from 12000 by 756 log-likelihood matrix.
## 
##          Estimate   SE
## elpd_loo   -697.7 32.0
## p_loo       157.8 10.7
## looic      1395.3 64.0
## ------
## MCSE of elpd_loo is NA.
## MCSE and ESS estimates assume MCMC draws (r_eff in [0.0, 0.9]).
## 
## Pareto k diagnostic values:
##                          Count Pct.    Min. ESS
## (-Inf, 0.7]   (good)     742   98.1%   14      
##    (0.7, 1]   (bad)       13    1.7%   <NA>    
##    (1, Inf)   (very bad)   1    0.1%   <NA>    
## See help('pareto-k-diagnostic') for details.
summarize_brms(
  reactance_ordinal, 
  model_rows_fixed = model_rows_fixed_ordinal,
  model_rows_random = model_rows_random_ordinal,
  model_rownames_fixed = model_rownames_fixed_ordinal,
  model_rownames_random = model_rownames_random_ordinal,
  exponentiate = T) %>%
  print_df(rows_to_pack = rows_to_pack_ordinal)
OR l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
Intercepts
Intercept NA NA NA NA NA NA
Intercept[1] 2.86* 1.41 6.43 1.005 1186.35 1593.65
Intercept[2] 6.64* 3.11 18.80 1.012 428.10 627.68
Intercept[3] 19.79* 8.32 81.74 1.017 299.84 325.99
Intercept[4] 96.91* 33.34 796.37 1.021 247.48 97.46
Intercept[5] 4844.88* 815.17 271440.98 1.023 228.76 112.06
Within-Person Effects
Daily persuasion experienced 1.03 0.90 1.20 1.001 5299.30 5209.39
Daily persuasion utilized (partner’s view) 0.91 0.75 1.08 1.003 3841.84 2192.73
Daily pressure experienced 0.87 0.66 1.14 1.003 4165.97 2964.96
Daily pressure utilized (partner’s view) 0.88 0.61 1.24 1.002 4402.75 4075.56
Daily pushing experienced 0.99 0.82 1.18 1.003 3554.27 2875.14
Daily pushing utilized (partner’s view) 1.08 0.87 1.34 1.003 4423.62 2792.88
Day 0.48 0.14 1.08 1.010 431.69 138.59
Daily weartime NA NA NA NA NA NA
Between-Person Effects
Mean persuasion experienced 0.88 0.32 2.29 1.006 942.94 446.93
Mean persuasion utilized (partner’s view) 1.12 0.38 3.29 1.005 905.01 884.69
Mean pressure experienced 2.01 0.77 5.47 1.002 4269.44 2934.53
Mean pressure utilized (partner’s view) 0.67 0.22 1.84 1.002 2613.68 2593.82
Mean pushing experienced 0.47 0.10 2.89 1.014 440.85 133.99
Mean pushing utilized (partner’s view) 1.26 0.26 8.85 1.011 441.83 131.88
Mean weartime NA NA NA NA NA NA
Random Effects
sd(Intercept) 1.34 0.91 2.16 1.01 288.61 155.72
sd(Daily persuasion experienced) 0.23 0.01 0.56 1.00 1682.37 2694.00
sd(Daily persuasion utilized (partner’s view)) 0.19 0.01 0.50 1.00 3004.62 5238.54
sd(Daily pressure experienced) 0.98 0.54 1.65 1.01 541.13 659.14
sd(Daily pressure utilized (partner’s view)) 0.41 0.02 1.22 1.00 2578.00 5800.65
sd(Daily pushing experienced) 0.21 0.01 0.60 1.00 1227.38 3317.84
sd(Daily pushing utilized (partner’s view)) 0.16 0.01 0.56 1.00 3278.27 5679.85
Additional Parameters
ar[1] 0.00 -0.14 0.14 1.00 12273.68 8802.84
nu NA NA NA NA NA NA
shape NA NA NA NA NA NA
sderr 0.78 0.02 2.84 1.03 164.15 102.07
sigma NA NA NA NA NA NA
disc 1.00 1.00 1.00 NA NA NA 
hypothesis(reactance_ordinal, "pressure_self_cw > pushing_self_cw")
## Hypothesis Tests for class b:
##                 Hypothesis Estimate Est.Error CI.Lower CI.Upper Evid.Ratio
## 1 (pressure_self_cw... > 0    -0.12      0.18    -0.41     0.17       0.32
##   Post.Prob Star
## 1      0.24     
## ---
## 'CI': 90%-CI for one-sided and 95%-CI for two-sided hypotheses.
## '*': For one-sided hypotheses, the posterior probability exceeds 95%;
## for two-sided hypotheses, the value tested against lies outside the 95%-CI.
## Posterior probabilities of point hypotheses assume equal prior probabilities.

Binary

introduce_binary_reactance <- function(data) {
  data$is_reactance <- factor(data$reactance > 0, levels = c(FALSE, TRUE), labels = c(0, 1))
  return(data)
}



df_double <- introduce_binary_reactance(df_double)
if (use_mi) {
  for (i in seq_along(implist)) {
    implist[[i]] <- introduce_binary_reactance(implist[[i]])
  }
}


formula <- bf(
  is_reactance ~ 
    persuasion_self_cw + persuasion_partner_cw +
    pressure_self_cw + pressure_partner_cw +
    pushing_self_cw + pushing_partner_cw +
    
    persuasion_self_cb + persuasion_partner_cb +
    pressure_self_cb + pressure_partner_cb +
    pushing_self_cb + pushing_partner_cb +
    
    day +
    
    # Random effects
    (persuasion_self_cw + persuasion_partner_cw +
    pressure_self_cw + pressure_partner_cw +
    pushing_self_cw + pushing_partner_cw | coupleID),

  autocor = ~ ar(time = day, gr = coupleID:userID, p = 1)
)



prior1 <- c(
  brms::set_prior("normal(0, 2.5)", class = "b"),
  brms::set_prior("normal(0, 10)", class = "Intercept", lb=0, ub=5), # range of the outcome scale
  
  brms::set_prior("normal(0, 2)", class = "sd", group = "coupleID", lb = 0),

  brms::set_prior("cauchy(0, 5)", class = "ar", lb = -1, ub = 1)
  #brms::set_prior("cauchy(0, 10)", class = "sigma", lb = 0)
)


#df_minimal <- df_double[, c("AorB", all.vars(as.formula(formula)))]

is_reactance <- my_brm(
  mi = use_mi, 
  imputed_data = implist,
  
  formula = formula, 
  prior = prior1,
  data = df_double, 
  family = bernoulli(),
  #control = list(adapt_delta = 0.95, max_treedepth = 15),
  iter = 5000,
  warmup = 2000,
  chains = 4,
  cores = 4,
  seed = 7777,
  file = file.path("models_cache_brms", "is_reactance")
)
plot(is_reactance, ask = FALSE)

plot(pp_check(is_reactance, type = 'ecdf_overlay'))
## Using 10 posterior draws for ppc type 'ecdf_overlay' by default.

plot(pp_check(is_reactance))
## Using 10 posterior draws for ppc type 'dens_overlay' by default.

#pp_check_transformed(is_reactance, transform = log1p)

loo(is_reactance)
## 
## Computed from 12000 by 756 log-likelihood matrix.
## 
##          Estimate   SE
## elpd_loo   -387.4 13.8
## p_loo       325.7 12.2
## looic       774.9 27.7
## ------
## MCSE of elpd_loo is NA.
## MCSE and ESS estimates assume MCMC draws (r_eff in [0.7, 1.2]).
## 
## Pareto k diagnostic values:
##                          Count Pct.    Min. ESS
## (-Inf, 0.7]   (good)       8    1.1%   471     
##    (0.7, 1]   (bad)      448   59.3%   <NA>    
##    (1, Inf)   (very bad) 300   39.7%   <NA>    
## See help('pareto-k-diagnostic') for details.
summarize_brms(
  is_reactance, 
  model_rows_fixed = model_rows_fixed,
  model_rows_random = model_rows_random,
  model_rownames_fixed = model_rownames_fixed,
  model_rownames_random = model_rownames_random,
  exponentiate = T) %>%
  print_df(rows_to_pack = rows_to_pack)
OR l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
Intercept 2.17 0.37 13.35 1.000 6516.79 7916.87
Within-Person Effects
Daily persuasion experienced 1.18 0.63 2.23 1.000 5629.41 7273.16
Daily persuasion utilized (partner’s view) 0.98 0.45 2.09 1.000 5509.28 6402.62
Daily pressure experienced 0.83 0.25 2.52 1.000 5560.23 6150.17
Daily pressure utilized (partner’s view) 0.95 0.21 4.23 1.001 6527.92 7616.08
Daily pushing experienced 0.99 0.46 2.08 1.000 5716.27 6967.70
Daily pushing utilized (partner’s view) 0.98 0.37 2.54 1.001 5715.73 6643.78
Day 0.18 0.01 2.98 1.000 7609.66 8520.45
Daily weartime NA NA NA NA NA NA
Between-Person Effects
Mean persuasion experienced 1.10 0.08 16.32 1.000 6229.05 7781.95
Mean persuasion utilized (partner’s view) 1.54 0.10 26.99 1.001 7151.80 8127.81
Mean pressure experienced 1.83 0.10 33.22 1.000 6720.79 7348.79
Mean pressure utilized (partner’s view) 0.89 0.05 16.38 1.000 7256.38 8303.77
Mean pushing experienced 0.13 0.00 5.00 1.000 7906.31 8605.17
Mean pushing utilized (partner’s view) 2.42 0.05 106.85 1.000 8609.03 8905.87
Mean weartime NA NA NA NA NA NA
Random Effects
sd(Intercept) 6.08 3.96 8.59 1.00 3875.52 5218.41
sd(Daily persuasion experienced) 1.19 0.17 2.48 1.00 1665.48 1795.99
sd(Daily persuasion utilized (partner’s view)) 1.14 0.07 2.61 1.00 1959.22 3003.34
sd(Daily pressure experienced) 3.14 1.47 5.31 1.00 4134.11 5651.54
sd(Daily pressure utilized (partner’s view)) 1.20 0.04 3.57 1.00 3961.08 5157.92
sd(Daily pushing experienced) 0.81 0.04 2.12 1.00 2291.64 4622.12
sd(Daily pushing utilized (partner’s view)) 0.74 0.03 2.22 1.00 4079.24 5495.08
Additional Parameters
ar[1] 0.11 -0.09 0.31 1.00 2089.64 4212.41
nu NA NA NA NA NA NA
shape NA NA NA NA NA NA
sderr 7.04 4.13 10.83 1.00 2189.05 3343.39
sigma NA NA NA NA NA NA 
hypothesis(is_reactance, "pressure_self_cw > pushing_self_cw")
## Hypothesis Tests for class b:
##                 Hypothesis Estimate Est.Error CI.Lower CI.Upper Evid.Ratio
## 1 (pressure_self_cw... > 0    -0.17      0.74    -1.43     1.02       0.71
##   Post.Prob Star
## 1      0.41     
## ---
## 'CI': 90%-CI for one-sided and 95%-CI for two-sided hypotheses.
## '*': For one-sided hypotheses, the posterior probability exceeds 95%;
## for two-sided hypotheses, the value tested against lies outside the 95%-CI.
## Posterior probabilities of point hypotheses assume equal prior probabilities.

Report All Models

if (report_ordinal) {
  model_rows_random_final <- model_rows_random_ordinal
  model_rows_fixed_final <- model_rows_fixed_ordinal
  model_rownames_fixed_final <- model_rownames_fixed_ordinal
  model_rownames_random_final <- model_rownames_random_ordinal
  rows_to_pack_final <- rows_to_pack_ordinal
} else {
  model_rows_random_final <- model_rows_random
  model_rows_fixed_final <- model_rows_fixed
  model_rownames_fixed_final <- model_rownames_fixed
  model_rownames_random_final <- model_rownames_random
  rows_to_pack_final <- rows_to_pack
}



all_models <- report_side_by_side(
  pa_sub,
  pa_obj_log,
  mood,
  reactance,
  is_reactance,
  
  model_rows_random = model_rows_random_final,
  model_rows_fixed = model_rows_fixed_final,
  model_rownames_random = model_rownames_random_final,
  model_rownames_fixed = model_rownames_fixed_final
) 
## [1] "pa_sub"
## [1] "pa_obj_log"
## [1] "mood"
## [1] "reactance"
## [1] "is_reactance"
# pretty printing
summary_all_models <- all_models %>%
  print_df(rows_to_pack = rows_to_pack_final) %>%
  add_header_above(
    c(" ", "Subjective MVPA" = 2, 
      "Device-Based MVPA" = 2, 
      "Mood" = 2,
      "Reactance Gaussian" = 2, 
      "Reactance Dichotome" = 2)
  )

export_xlsx(
  summary_all_models, 
  rows_to_pack = rows_to_pack_final,
  file.path("Output", "AllModels.xlsx"), 
  merge_option = 'both', 
  simplify_2nd_row = TRUE,
  colwidths = c(38, 7.2, 13.3, 7.2, 13.3,7.2, 13.3,7.2, 13.3,7.2, 13.3),
  line_above_rows = c(1,2),
  line_below_rows = c(-1)
)

summary_all_models
Subjective MVPA
Device-Based MVPA
Mood
Reactance Gaussian
Reactance Dichotome
IRR pa_sub 95% CI pa_sub exp(Est.) pa_obj_log 95% CI pa_obj_log b mood 95% CI mood b reactance 95% CI reactance OR is_reactance 95% CI is_reactance
Intercept 24.25* [17.04, 34.72] 116.98* [104.97, 130.74] 4.93* [ 4.71, 5.15] 0.64* [ 0.37, 0.90] 2.17 [0.37, 13.35]
Within-Person Effects
Daily persuasion experienced 0.97 [ 0.87, 1.07] 1.00 [ 0.98, 1.02] 0.00 [-0.03, 0.03] 0.01 [-0.04, 0.06] 1.18 [0.63, 2.23]
Daily persuasion utilized (partner’s view) 1.05 [ 0.95, 1.17] 1.00 [ 0.98, 1.02] -0.01 [-0.05, 0.02] -0.04 [-0.11, 0.02] 0.98 [0.45, 2.09]
Daily pressure experienced 1.12 [ 0.89, 1.45] 1.03 [ 0.98, 1.09] 0.00 [-0.08, 0.07] -0.04 [-0.13, 0.06] 0.83 [0.25, 2.52]
Daily pressure utilized (partner’s view) 0.99 [ 0.79, 1.28] 0.98 [ 0.93, 1.03] 0.05 [-0.03, 0.13] -0.04 [-0.17, 0.09] 0.95 [0.21, 4.23]
Daily pushing experienced 0.93 [ 0.80, 1.09] 1.00 [ 0.97, 1.03] -0.03 [-0.08, 0.02] -0.01 [-0.07, 0.05] 0.99 [0.46, 2.08]
Daily pushing utilized (partner’s view) 0.92 [ 0.79, 1.09] 1.02 [ 0.99, 1.06] -0.04 [-0.09, 0.00] 0.02 [-0.06, 0.10] 0.98 [0.37, 2.54]
Day 1.30 [ 0.95, 1.78] 1.03 [ 0.94, 1.12] -0.10 [-0.25, 0.06] -0.25 [-0.52, 0.03] 0.18 [0.01, 2.98]
Daily weartime NA NA 1.00 [ 1.00, 1.00] NA NA NA NA NA NA
Between-Person Effects
Mean persuasion experienced 1.21 [ 0.72, 2.06] 1.05 [ 0.90, 1.22] 0.14 [-0.13, 0.42] -0.04 [-0.39, 0.31] 1.10 [0.08, 16.32]
Mean persuasion utilized (partner’s view) 1.16 [ 0.71, 1.94] 1.05 [ 0.90, 1.23] 0.21 [-0.06, 0.49] 0.09 [-0.30, 0.47] 1.54 [0.10, 26.99]
Mean pressure experienced 1.60 [ 0.82, 3.12] 0.86 [ 0.72, 1.02] -0.11 [-0.37, 0.14] 0.29 [-0.07, 0.65] 1.83 [0.10, 33.22]
Mean pressure utilized (partner’s view) 0.53 [ 0.28, 1.04] 0.87 [ 0.74, 1.01] -0.02 [-0.28, 0.23] -0.22 [-0.60, 0.16] 0.89 [0.05, 16.38]
Mean pushing experienced 0.43* [ 0.19, 0.99] 1.12 [ 0.90, 1.40] -0.20 [-0.60, 0.21] -0.24 [-0.77, 0.30] 0.13 [0.00, 5.00]
Mean pushing utilized (partner’s view) 0.68 [ 0.29, 1.55] 1.12 [ 0.90, 1.38] -0.43* [-0.83, -0.03] 0.03 [-0.56, 0.60] 2.42 [0.05, 106.85]
Mean weartime NA NA 1.00 [ 1.00, 1.00] NA NA NA NA NA NA
Random Effects
sd(Intercept) 0.68 [ 0.51, 0.91] 0.29 [0.22, 0.38] 0.63 [0.49, 0.81] 0.57 [ 0.43, 0.76] 6.08 [ 3.96, 8.59]
sd(Daily persuasion experienced) 0.22 [ 0.07, 0.37] 0.06 [0.03, 0.09] 0.03 [0.00, 0.07] 0.06 [ 0.00, 0.14] 1.19 [ 0.17, 2.48]
sd(Daily persuasion utilized (partner’s view)) 0.18 [ 0.03, 0.33] 0.05 [0.02, 0.08] 0.06 [0.01, 0.11] 0.04 [ 0.00, 0.11] 1.14 [ 0.07, 2.61]
sd(Daily pressure experienced) 0.14 [ 0.01, 0.44] 0.05 [0.00, 0.14] 0.09 [0.00, 0.26] 0.44 [ 0.28, 0.67] 3.14 [ 1.47, 5.31]
sd(Daily pressure utilized (partner’s view)) 0.15 [ 0.01, 0.45] 0.04 [0.00, 0.11] 0.15 [0.01, 0.33] 0.22 [ 0.02, 0.58] 1.20 [ 0.04, 3.57]
sd(Daily pushing experienced) 0.25 [ 0.02, 0.53] 0.08 [0.01, 0.15] 0.09 [0.02, 0.16] 0.09 [ 0.00, 0.23] 0.81 [ 0.04, 2.12]
sd(Daily pushing utilized (partner’s view)) 0.14 [ 0.01, 0.36] 0.04 [0.00, 0.10] 0.08 [0.01, 0.16] 0.04 [ 0.00, 0.14] 0.74 [ 0.03, 2.22]
Additional Parameters
ar[1] 0.02 [-0.93, 0.94] 0.28 [0.25, 0.32] 0.45 [0.42, 0.48] 0.01 [-0.07, 0.10] 0.11 [-0.09, 0.31]
nu NA NA NA NA NA NA NA NA NA NA
shape 0.14 [ 0.13, 0.14] NA NA NA NA NA NA NA NA
sderr 0.05 [ 0.00, 0.13] NA NA NA NA NA NA 7.04 [ 4.13, 10.83]
sigma NA NA 0.55 [0.54, 0.57] 0.87 [0.85, 0.89] 0.94 [ 0.89, 0.99] NA NA 
report::report_system()

Analyses were conducted using the R Statistical language (version 4.4.1; R Core Team, 2024) on Windows 11 x64 (build 22635)

report::cite_packages()